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We observe using ab initio methods that localized surface plasmon resonances in icosahedral silver 
nanoparticles enter the asymptotic region already between diameters of 1-2 nm, converging close to 
the classical quasistatic limit around 3.4 eV. We base the observation on time-dependent density- 
functional theory simulations of the icosahedral silver clusters Agss (1.06 nm), Agi 47 (1.60 nm), 

Ag 309 (2.14 nm), and Agsgi (2.68 nm). The simulation method combines the adiabatic GLLB-SC 
exchange-correlation functional with real time propagation in an atomic orbital basis set using the 
projector augmented wave method. The method has been implemented to the electron structure 
code GPAW within the scope of this work. We obtain good agreement with experimental data and 
modelled results, including photoemission and plasmon resonance. Moreover we can extrapolate the 
ab initio results to the classical quasistatically modelled icosahedral clusters. 


I. INTRODUCTION 

Localized surface plasmon resonances (LSPR) of sil¬ 
ver nanoparticles (AgNPs) exhibit strong UV-VIS ab¬ 
sorption. The LSPR can be tuned by fabrication 
techniques,^ or by functionalization,^ and it is sensitive 
to the nanoparticle’s environment.^ Sensitivity and tun- 
ability of AgNPs can be utilized in sensing,^ surface- 
enhanced spectroscopies,® plasmon-enhanced chemistry,® 
and photovoltaic applications.^ Much of the wide inter¬ 
est in AgNPs originates from their role as building blocks 
of nanophotonic devices, such as optical nanoantennas.® 
The ability to predict the relation between their structure 
and operation is crucial for the applications. The optical 
characteristics of large noble metal NPs (> 10 nm) are 
well known, and their LSPR can be simulated using clas¬ 
sical electromagnetic theory. For example, large spherical 
AgNPs have a LSPR at 355 nm (3.5 eV), whereas icosa¬ 
hedral particles are slightly redshifted and have broader 
absorption arising from several LSPR modes that overlap 
closely in energy.^do However, as the diameter of the NPs 
decreases, the LSPR blueshifts with the frequency being 
inversely proportional to the diameter^ ^ and finally, the 
absorption spectrum changes to a typical cluster spec¬ 
trum characterized by several individual transitions be¬ 
tween quantized energy levels.For diameters smaller 
than 10 nm, the sensitivity of the LSPR to the shape and 
surroundings of the AgNP becomes important which is 
reflected in the difficulty of interpretation of experiments. 

Previous theoretical studies on the LSPR in AgNPs 
are limited to quantum mechanical calculations of small 
clusters^"*’ and jellium models,^® or to classical electro¬ 
magnetic theory for large NPs.® Between diameters of 
1-5 nm, the classical electromagnetic theory does not 


provide an adequate description of the possible quantum 
effects as it is scale invariant and therefore predicts no 
size dependence. Jellium models ignore — or at best 
approximate^® — the effect of d-electrons and atomic 
structure which are crucial to the proper description of 
AgNPs. Ab initio methods are limited to small clusters: 
Time-dependent^^ (TD) density-functional theory^®d9 
(DFT) has been used to model Agss and Agi 47 ,^^ as well 
as nanoshells up to Ag 272 .^® 

Typically in such studies one uses the adiabatic lo¬ 
cal density approximation (ALDA) or adiabatic gen¬ 
eralized gradient approximation (AGGA) as exchange- 
correlation (XG) functionals, even though LDA and GGA 
are known to predict a too high-lying d-electron band and 
therefore to severely overestimate the d-band screening. 
This results in decreased oscillator strength and lowered 
plasmonic frequency compared to experiments. Both 
experimental^^ and theoretical^® works have conhrmed 
that the position of the d-band strongly influences the 
plasmonic properties. Quantitative theory must be based 
on a more accurate description of the d-band. 

A recent experimental study of Scholl et al. found 
quantum effects influencing the optical properties of Ag¬ 
NPs with diameters as large as 10 nm.®"^ In particular, 
their electron-energy-loss spectroscopic (EELS) measure¬ 
ments on NPs showed a significant (0.5 eV) blueshift of 
the LSPR when the diameter decreased from 7 nm to 
2 nm. This disagrees with previous experimental results 
for freestanding clusters,®®’^® and Haberland has sug¬ 
gested that the blueshift is not due to the quantum effects 
but due either to the interaction of the LSPR with the 
substrate or the residual ligand molecules.®® This con¬ 
troversy exemplifies that without tools that can simulate 
the optical properties of NPs from molecular size up to 
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the classical limit, it is difhcult to separate the quantum 
effects from other factors. 

In this work, we present an ab initio theoretical anal¬ 
ysis of freestanding AgNPs up to diameter of 2 nm, and 
show that it is unlikely that the results of Scholl et al. 
correspond to freestanding AgNPs. Using atomistic first- 
principles calculations in the TDDFT framework, we are 
able to obtain the macroscopic LSPR already at a di¬ 
ameter of 2 nm. Our calculations show that the res¬ 
onance shifts only by 0.2 eV above that. Our results 
agree with the experimental cluster data for both the 
smallest and the largest structures. Concurrently with 
explaining the experimental findings, we show that ac¬ 
curate treatment of interband (d-electron) excitations is 
crucial for a reliable description of AgNP plasmonics. 
Therefore, we recommend the adiabatic Gritsenko-van 
Leeuwen-van Lenthe-Baerends—solid-correlation poten¬ 
tial (GLLB-SC)^^ for approximating the exchange and 
correlation effects for the optical properties of noble 
metal NPs. The potential is a modification of the GLLB- 
potentiaU® to be better suited for solids and surfaces and 
with added correlation. 

In Section II we describe the details of our implementa¬ 
tion of linear combinations of atomic orbitals with time- 
dependent density-functional theory (LCAO-TDDFT), 
and elaborate the relevance of the GLLB-SG potential 
for the proper description of plasmonics. In Section III 
we give basic background information about the quan¬ 
tum mechanical and the electrodynamical model. In Sec¬ 
tion IV we analyse the obtained results, and compare 
them to experimental EELS and photoemission data. In 
Section V we carefully benchmark the accuracy of our 
method. In Section VI we summarize the results and 
discuss the relevance of proper Kohn-Sham eigenvalue 
description for accurate absorption spectra in AgNPs. 

II. METHODS 

The main computational challenges in simulating the 
photoabsorption spectrum of nanoplasmonic structures 
using TDDFT are 1) quality of the XC functional, espe¬ 
cially for the description of the silver d-band, 2) the nu¬ 
merical discretization scheme for the wavefunctions and 
the density, which must be flexible enough to describe the 
LSPR, and 3) the method for optical properties must be 
fast, parallelizable, and scale well with respect to sys¬ 
tem size. Each of the challenges will be addressed in the 
following subsections. 

A. Time-dependent Density-Functional Theory 

Time-dependent density-functional theory is a well es¬ 
tablished tool for calculating electronic excitations. As in 
DFT, the most crucial aspect of TDDFT is the exchange- 
correlation potential, which is time-dependent in this 
case. The time-dependent Kohn-Sham equations for the 


electronic orbitals 'I'i are 

^-i9t - -b■^;Ks[?^(l'U)](l■U)^ ^i(rU) = 0, (1) 

where uks is the Kohn-Sham potential, the time- 
dependent density is given by 

n{r,t) (2) 

i 

and fi are the occupation numbers of the orbitals. 

In the general formalism, the exchange-correlation 
part of the Kohn-Sham potential t>xc depends causally 
on all previous densities. In a practical and widely used 
adiabatic approximation, the potential depends only on 
the instantaneous density. We will use this approxima¬ 
tion also in the case of the GLLB-SG potential, with one 
further modification, as discussed in the next section. 


B. Adiabatic GLLB-SC 

Adiabatic (semi)local density approximations, such as 
ALDA and AGGA, are applicable for nearly free-electron 
metals, but for noble metals the situation is differ¬ 
ent because they overestimate the polarizability of d- 
electrons. This is due to their Kohn-Sham spectrum, 
since they predict too delocalized d-band in addition 
it being too shallow.To overcome this problem we 
employ the adiabatic GLLB-SC potentiaP^’^® that in¬ 
cludes the exchange-hole and correlation potential of the 
Perdew-Burke-Ernzerhof functional for solids and sur¬ 
faces (PBEsol),®° and is additionally supplemented by a 
computationally efficient approximation of the hole re¬ 
sponse part (see, e.g.. Ref. 31) of the exact-exchange op¬ 
timized effective potential.®^ 

GLLB-SC introduces an orbital energy dependent 
localization of the exchange hole which reduces self¬ 
interaction and yields better asymptotic behavior than 
LDA or GGA.^^ So far GLLB-SC has been mostly ap¬ 
plied for predicting semiconductor band gaps,^^’®® but 
recently Yan et al. showed that it also yields good re¬ 
sults for Ag surface plasmons because of the improved 
d-band description.^^ We employ this finding, but ex¬ 
tend it further by applying GLLB-SC also for the dy¬ 
namic response (in Ref. 21, GLLB-SC was used only for 
the ground state whereas the linear response calculation 
employed ALDA). 

We obtain an adiabatic GLLB-SC approximation 
by replacing the time-dependent response coefficients 
Wi{t) with their time-independent ground-state values 
Wi{t = 0) = Kgy/€f — Ci in the GLLB-SC potential (see 
Eqns. (16) and (22) of Ref. 27): 

I'GLLB-SC(I’) = Vx-hole(r) + y^ Wi{t = 0) -b 'Uc(r), 

nr) 

i 

( 3 ) 
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where Ux-hoie(r) is the Coulomb potential due to the ex¬ 
change hole obtained from the exchange hole of PBEsol 
evaluated at the instantaneous density, Uc(r) is the semi¬ 
local PBEsol correlation potential, and the remaining 
term is an approximation to the response of the Coulomb 
potential of the exchange-correlation hole to density per¬ 
turbations. We choose Wi{t) = Wi{t = 0), since it is the 
simplest obtainable approximation and computationally 
attractive. It is plausible that this approximation is ac¬ 
curate in our simulations because we apply a small per¬ 
turbation which will not significantly change the density, 
and thus not induce large oscillations of Wiit). In addi¬ 
tion, our preliminary adiabatic time-dependent Krieger- 
Li-Iafrate (TD-KLI)^'^ calculations indicate that the ef¬ 
fect of Wiit) compared to Wiit = 0) in the systems con¬ 
sidered here are negligible in the linear response regime. 

It is difficult to estimate the effect of this approxima¬ 
tion, or the effect of non-adiabatic exchange-correlation 
effects exactly. However, there is much evidence from 
historical work that already the random phase approx¬ 
imation (pure Coulomb kernel) without any XC kernel 
is sufficient to describe plasmonics.^^ Therefore, the adi¬ 
abatic GLLB-SC approximation should be a sufficient 
description for AgNP plasmonics. 

C. Real Time Propagation with Basis Sets 

The wavefunctions are represented as linear combinations 
of atomic orbitals (LCAO) together with the projector 
augmented wave method^® (PAW) as implemented in the 
GPAW package.The smooth pseudo wavefunctions 
are written as a linear combination 


operator T are constant because the nuclei are assumed 
to be stationary. 

In this approach, the time-dependent density and po¬ 
tential are expressed on a uniform grid, and the matrix el¬ 
ements of the potential are evaluated on this grid.^® The 
smoothness of these quantities allows for a very coarse 
grid spacing, and the LCAO-PAW pseudo wavefunctions 
form a small, local, and efficient representation suitable 
for systems with hundreds of atoms.^® 

We calculate the optical absorption spectrum of Ag- 
NPs using the time-propagation (TP) approach to 
TDDFT."^*^di The greatest advantage of TP-TDDFT is 
the scaling of the computational requirements with re¬ 
spect to system size compared to other methods, such 
as Gasida’s approach.Despite its better scaling, the 
large prefactor has so far limited the applicability of the 
TP-TDDFT approach. 

Following the TP-TDDFT procedure for the optical 
response,we here excite the system by an instantaneous 
electric field E(r,t) = EoSttSit), where the field strength 
Eq — 0.0001 a.u. is sufficiently small to avoid nonlinear 
effects, and the direction ett of the electric field is chosen 
to be from tip to tip, i.e., along the five-fold symmetry 
axis of the icosahedron. The optical absorption spectrum 
is obtained by Fourier transforming the induced dipole 
moment along the excitation axis.^*^ 

After the initial kick, the propagation is performed 
with a reliable and numerically stable semi-implicit 
Grank-Nicolson method. In brief, the method can be 
described as follows. In the prediction step, we solve 

(^S + if H(t)^ ait + dt) =(s- if H(t)^ C(t), (7) 


^,(r,t) = ^C^,(t)0^(r-R^) (4) 

of atom centered orbitals ipnir — R^) with expansion 
coefficients C'^i(t). The PAW projection operator^® T 
can be used to reconstruct the all-electron functions 
as 4'i(r,t) = T^i(r,t). The PAW form of the time- 
dependent Kohn-Sham equations (1) is 


for C'(f-l-dt), where S and H(t) are the basis set represen¬ 
tations of the PAW overlap and Hamiltonian operator re¬ 
spectively. The operations are parallel with matrices be¬ 
ing distributed using ScaLAPAGK^^ and BLAGS.^® Af¬ 
ter obtaining the initial approximation for the wavefunc¬ 
tions, the predict-correct method is applied. We then 
obtain an estimate for the Kohn-Sham Hamiltonian (in¬ 
cluding the XG potential) at the middle of the time step 


TM-i-)f+ f^RKs(t)f 


T,(r,t) = 0, (5) 


where Rks( 0 is the Kohn-Sham Hamiltonian for non¬ 
interacting electrons. 

Substituting Eq. (4) into Eq. (5) and multiplying with 
J dr(^^(r) from the left, the equation can be cast into a 
matrix form 

iS^=H(f)C(f), (6) 

with the overlap matrix 5^,^ = ((^^jT'TIt^jy) and the 
Hamiltonian matrix H^i,it) = H]^s,it)T\^u)■ C(t) 

is the matrix of LG AO expansion coefficients {C'^i(t)} 
defined in Eq. (4). The overlaps and the projection 


H(t + dV2)«(H(f) + H'(t + df))/2, (8) 

where H'(t -f dt) is evaluated from C'(t -|- dt), and then 
propagate the wavefunction to t -I- dt in the correction 
step which solves 

-f if H(t -b dt/2)^ C(t -f dt) 

= (^S-ifH(t + dt/2)^ C(t), (9) 

for C(t-l-dt). This results in C>(A®) scaling with respect 
to the number of electrons in the system, to be compared 
to the GPAW’s Gasida implementation of OiN^) or the 
real-space time propagations C>(7V^). However, the con¬ 
stant factor in the grid propagation is so large, that our 
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scheme performs 1 to 3 orders of magnitude faster on 
systems of several thousand electrons. The timings for 
propagation are indicated in Table I. 


TABLE I. Performance of the LCAO-TDDFT code for time 
propagation of AgNPs for 1000 time steps of duration 10 as. 


System 

Cores Wall hrs CPU hrs Electrons Basis functions 

Agss 

64 

4.5 

288 

605 

990 

Agl47 

64 

18.0 

1152 

1617 

2646 

Ag309 

256 

28.5 

7296 

3399 

5562 

Ag561 

512 

42.0 

21504 

6171 

10098 


III. MODEL 

Both classical and quantum mechanical models are em¬ 
ployed in this work. The quantum mechanical model is 
atomistic and ab initio, relying only on DFT and TDDFT 
calculations. The plasmonic peak in our quantum me¬ 
chanical model depends both on shape and size effects. 
The classical model is based on empirical dielectric func¬ 
tions and can only model shape effects. However, in 
the large particle limit, the two methods should agree. 
Therefore, we can test the performance of our computa¬ 
tional model also by extrapolating to macroscopic Mie 
scattering limit. 

In both models we consider icosahedral clusters. 
Charged Agss has been experimentally identified as 
icosahedral,^® and we choose to keep the icosahe¬ 
dral geometry to avoid shape effects even though the 
minimum energy structure is expected to change for 
larger clusters.The difference between icosahedron and 
sphere in the macroscopic limit is well understood®hO and 
does not influence our conclusions. The atomic structure 
of icosahedral Agss used in our calculations includes a 
central atom with two icosahedral Mackay layers. We 
create clusters up to Agsgi by adding Mackay layers one 
by one, using a bond length of 3.0 A. The ideal icosahe¬ 
dral clusters are then relaxed with the LDA functional. 
In these calculations the grid spacing is 0.2 A, and the 
size of the cubic cell is chosen so that all atoms are at 
least 5.0 A away from the cell boundary. We use the de¬ 
fault double-C polarized (DZP) basis set provided with 
GPAW for the geometry relaxations.®® 

TDDFT simulations are performed for 30 fs using time 
steps of 10 as. These calculations use a coarse grid spac¬ 
ing of 0.3 A and an expanded atomic basis set; we will 
further discuss these parameters in Section V. All spec¬ 
tra are calculated using a Gaussian broadening of 0.16 eV 
FWHM. 

In the following we consider classical electrodynamic 
approximations. First, the photoabsorption of a spherical 
NP of volume V is given by the quasistatic limit of Mie 


theory: 


By using the experimentally determined permittivity 
e(a;) for silver presented in Ref. 48, Eq. (10) yields a 
strong LSPR at 3.5 eV. For more complicated shapes, 
such as icosahedra, one has to employ computational 
electrodynamics. In this work we use a quasistatic 
(QS) version^® of the widely used finite-difference-time- 
domain (FDTD) method, as implemented in GPAW.®® 
Like in photoabsorption calculations with TDDFT, in 
the QSFDTD method one perturbs the system by an 
external held and analyses the time-dependent dipole 
moment. The frequency-dependent dielectric permittiv¬ 
ity of classical material is approximated using a set of 
Lorentzians. To obtain an accurate representation of the 
dielectric function of Ag especially near the LSPR, we 
start from the parametrization presented in Ref. 49 which 
uses 9 Lorentzians, add one extra Lorentzian, and reht 
the dielectric function against the experimental data^® 
with weight function 10 ( 10 ) = exp(—(w/eV — 3.5)^). 

The QSFDTD calculations are performed using a reg¬ 
ular grid of 96 X 96 X 96 points. Since this method is 
size invariant, only a particle shape needs to be speci- 
hed. We thus specify the shape as an icosahedron with 
a length of 40 points along its axis, securing adequate 
surrounding vacuum. The material is represented by a 
mask which assigns a value of either 1 (material) or 0 
(vacuum) to each point. To ensure high numerical accu¬ 
racy of the finite-difference operators, we smoothen the 
edge of the icosahedron artificially over 2-3 grid points 
along the faces so that points along the faces are effec¬ 
tively a mixture of vacuum and silver. 

IV. RESULTS 

Fig. 1(a) shows the GLLB-SG TDDFT absorption spec¬ 
tra of icosahedral Agss, Agi 47 , Agsog, and Agsei clusters 
divided by number of atoms in the system. For com¬ 
parison, we present classical QSFDTD results for icosa¬ 
hedral (dashed red) and spherical shape (dashed black). 
These correspond to the limit of large clusters as given by 
the quasistatic approximation. Fig. 1(b) shows excitation 
energies of absorption peaks with respect to the inverse 
diameter of the cluster. For NPs larger than Ag 55 , the 
excitation energies of the most intense peak as a function 
of inverse diameter lie on a line (solid red) which extrap¬ 
olates to 3.35 eV in the large particle limit, very close 
to the mesoscopic limit for icosahedral AgNPs at 3.43 eV 
(dashed red) obtained from the QSFDTD calculation. 
The agreement of the quasistatic mesoscopic limit with 
the quantum mechanical asymptotic limit suggests that 
the quantum mechanical model correctly describes the 
shape effect. For comparison, a linear fit to experimental 
data (solid black) is shown for spherical AgNPs in argon 
matrix®^ and also the mesoscopic limit for spherical NPs 


Siu:) = ^Im 
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FIG. 1. a) Photoabsorption spectra for AgNPs normalized by 
number of atoms. The spectra for icosahedral Agss, Agi 47 , 
Ag 309 and Agsei are calculated with adiabatic GLLB-SC 
TDDFT. For comparison, the PBE-calculated spectrum of 
Agi 47 is also shown. Glassical QSFDTD spectra are shown 
for spherical (QS Sph) and icosahedral (QS Ico) AgNPs, cal¬ 
culated with the empirical dielectric function. These spectra 
are normalized using the empirical silver density and thus 
the peak strengths are directly comparable, b) LSPR energy 
(circles) of icosahedral AgNPs as a function of inverse par¬ 
ticle diameter with GLLB-SG TDDFT. Solid black and red 
lines correspond to the experimental data for spherical NPs^^ 
and a linear fit to our results, respectively. Dashed horizontal 
lines represent classical limits (QSFDTD) for spherical (QS 
Sph) (3.52 eV) and icosahedral (QS Ico) (3.43 eV) NPs. 


from the QSFDTD calculation (dashed black). The ex¬ 
perimental values are shifted to the vacuum LSPR value 
of 3.5 eV to account for the Ar matrix as suggested by 
Haberland.^® The experimental and the simulated data 
show remarkable agreement, both in the asymptotic limit 
and in the size dispersion. The differences can be at¬ 
tributed to slightly different AgNP shape and structure. 
These observations suggest that the quantum mechan¬ 
ical model describes the finite size effect in the AgNP 
plasmonics well. 

In Fig. 1(a), in addition to the LSPR energy, also the 
area of the plasmon peak per particle (oscillator strength) 
agrees well with the classical electrodynamics simulation 
(dashed red). These observations strongly indicate that 
(i) adiabatic GLLB-SC provides realistic d-band screen¬ 
ing in Ag nanostructures, and (ii) the macroscopic size 
range is reached for AgNPs of diameter ^2 nm. In Fig. 1, 
for comparison, we have included the spectrum of Agi 47 
calculated with PBE.®^ Importantly, the previous conclu¬ 
sions cannot be drawn from APBE calculations because 
they underestimate the LSPR energy by ~0.5 eV, and 
greatly underestimate the intensity as seen on Fig. 1(a). 

Previous works®^’^^ have demonstrated the importance 
of visual interpretation for characterizing the LSPR in 
molecules and NPs. The induced electron densities of 
LSPRs in Ag 55 , Agi 47 , Agaog, and Agsei are shown in 


Fig. 2. The exact quantity shown is the transition density 
at the plasmon frequency w of each AgNP, i.e., a sine 
transform 

nOO 

fi{r,uj) = / dt [n(r, t) — n(r, 0)] e“°’* sinwt (11) 

Jo 

of the charge density fluctuation. The damping is given 
by ctfwhm = 2-^210^(7 = 0.16 eV. We observe that 
the Ag sp-band near the Fermi energy forms a localized 
surface plasmon mainly at the two opposing sides of the 
icosahedron, whereas d-electrons polarize in the opposite 
direction and thus create a counteracting screening field 
at the central region. This screening is overestimated by 
PBE, causing the drop in plasmon energy and intensity. 
The figure corresponds to the classical picture of plas- 
mons as a charge cloud oscillating between the opposite 
sides of the AgNP. The visual inspection thus supports 
our finding that the macroscopic plasmon forms in the 
clusters of this size range. 



FIG. 2. Galculated induced electron densities of LSPR for 
a) Ag 55 , b) Agi47, c) Ag309, and d) Agsei. The Ag sp-band 
forms a localized surface plasmon on the surface of the cluster, 
whereas the d-electrons polarize in the opposite direction. 

Fig. 3 shows the experimental photoemission data from 
two sources®^’®® on Agss compared to sp and d-band pro¬ 
jected local density of states of the quantum mechanical 
clusters. The d-band position of GLLB-SC matches well 
with the experimental data, as has also been observed 
earlier.In addition, the superatom shell description is 
in quantitative agreement with photoemission data.^^ 


V. ACCURACY OF THE METHOD 

The PAW dataset used to represent Ag includes the 5s 
and 4d orbitals as valence states, and is based on the 
default parameters of GPAW for the 11-electron Ag setup 
(e.g., the PBE Ag setup from GPAW-setups vO.8.7929) 
but generated with the GLLB-SC functional. 
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FIG. 3. Density of states of icosahedral AgNPs calculated 
with GLLB-SC, projected separately for d-band (red) and 
sp-band (blue). The black curves are experimental data,®"*’®® 
shifted to align the Fermi levels. The charge state of Agss has 
insignificant effect on the quantitative agreement. 


In GPAW, one commonly uses a DZP numerical basis 
set to represent the wavefunctions.^® This basis set in¬ 
cludes the atomic Kohn-Sham orbital for each occupied 
valence state, one extra radial function for each atomic 
KS orbital generated using the standard “split-valence” 
scheme in GPAW, plus a polarization function which for 
transition metals is p-type. For the details on the con¬ 
struction of the basis sets, see Ref. 38. This basis set 
is designed for ground-state calculations and would not 
be expected to (and indeed does not) accurately predict 
properties that depend of unoccupied states. To better 
represent the effect of the unoccupied 5p orbitals, we re¬ 
place the standard p-type polarization function with the 
actual Kohn-Sham orbital of the 5p state plus its usual 
split-valence function. 



Energy (eV) 

FIG. 4. The density of states of (a) the Agss cluster and (b) 
the Agi 47 duster calculated with different LGAO basis sets 
as well as grid mode. 


Outside of this, we use the specific generation 
parameters®® of 0.07 eV confinement energy to localize 
the KS orbitals and a tail norm of 0.2 to define the range 
of the split-valence functions. These latter parameters we 
have optimized to provide an accurate density of states 
(DOS) in Agss as compared to an accurate real-space 
grid calculation, but this optimization has very little ef¬ 
fect compared to the inclusion of the diffuse 5p valence 
orbital. A comparison of DOSs is presented in Fig. 4. 
We observe that without the diffuse 5p valence orbital 
the basis set is not able to reproduce the correct DOS 
accurately, particularly for high energies. 

Fig. 5 presents the photoabsorption spectrum of the 
Agi 47 cluster calculated with different basis sets and on 
a real-space grid. Like in the DOS comparison, we note 
that the enhanced basis yields significantly better agree¬ 
ment with the grid mode than the default basis sets. In 
comparison to the real-space calculation, the enhanced 
LCAO basis reproduces the spectrum within ^0.1 eV 
and ^ 5% accuracy for peak energy and intensity, re¬ 
spectively. This approach yields a transferable basis set 
that can be expected to describe both DOS and the op¬ 
tical response of larger clusters with good accuracy. To 
obtain further improvement in accuracy, more elaborate 
approaches can be used to enhance the basis set.®® 



Energy (eV) 


FIG. 5. The photoabsorption spectrum of the Agi 47 cluster 
calculated with different LCAO basis sets and the grid mode. 


To obtain good convergence with respect to the vac¬ 
uum size, it is essential not to use zero boundary con¬ 
ditions for solving the Hartree potential of the Pois¬ 
son equation. In the current work, we employ a mul¬ 
tipole moment expansion®^ in order to obtain the correct 
boundary values of the Hartree potential. This allows us 
to scale down the required amount of vacuum from 15 A 
to 5 A and obtain significant speedup. 

As indicated by Table I, our method achieves good 
parallel scaling in the weak sense, i.e., the computational 
time can be kept within reasonable limits by increasing 
the number of CPU cores as the system size increases. 
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VI. CONCLUSIONS 

We have established and carefully benchmarked a real 
time propagation method using atomic basis sets to ob¬ 
tain accurate plasmonics as demonstrated here for icosa- 
hedral silver clusters. The implementation is part of 
the free open-source GPAW package. We have shown 
that the eigenvalue spectrum of the GLLB-SC potential 
matches the available experimental photoemission data 
for icosahedral silver clusters and that the method pro¬ 
vides an accurate description of the plasmonic response 
in TDDFT calculations. 

The observation that only the LSPR of Agss does not 
fit the asymptotic line in Fig. 1(b) suggests that the 
macroscopic regime is reached already at Agi 47 . How¬ 
ever, comparison of the spectrum of Agi 47 with the larger 
clusters shows that the shape of the LSPR peak deviates 
from the larger clusters. These quantum effects disap¬ 
pear for Ag 309 and larger clusters. This threshold size 
for asymptotic LSPR behavior is remarkably small and 
agrees with experimental observations^^’^® as well as with 
simulations of monolayer protected Au clusters.®^ 

The impact of this study is threefold. Firstly, we show 
using ab initio simulations that the LSPR frequencies 
and intensities in icosahedral AgNPs enter an asymp¬ 
totic region already around the diameter of 2 nm. The 
optical response converges close to the classical limit of 
3.43 eV for icosahedral AgNPs. Our simulations are in 
good agreement with the experimental data, and the con¬ 
clusion is further supported by visual examination and 
analysis of the DOS. The presented results thus set the 
benchmark for the plasmonics of AgNPs, and explain 
the controversy between the recent EELS results with 


previous cluster experiments.Secondly, the results 
show that adiabatic GLLB-SC provides an accurate de¬ 
scription of d-band screening in Ag nanostructures with 
computational effort that is comparable to ALDA and 
AGGAs. The final point of the study — with probably 
the greatest impact in the long run — is the efficiency 
of the combination of TP-TDDFT, LCAO and the PAW 
method. The method is not limited to pure Ag nanos¬ 
tructures. Our preliminary results show that it is also ap¬ 
plicable to intermetallic nanostructures, such as Au~Ag 
core-shell NPs, as well as to nanostructures with molec¬ 
ular parts, e.g., ligand protected AuNPs®^ and metallic 
nanoantennas connected by molecular tunnel junctions.^ 
Altogether the combination of adiabatic GLLB- 
SC, LCAO-PAW with extended basis, and the time- 
propagation method allows for accurate simulations of 
LSPRs in noble metal nanostructures towards macro¬ 
scopic sizes. 
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